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Abstract Inversions for local helioseismology are an important and necessary 
step for obtaining three-dimensional maps of various physical quantities in the 
solar interior. Frequently, the full inverse problems that one would like to solve 
prove intractable because of computational constraints. Due to the enormous 
seismic data sets that already exist and those forthcoming, this is a problem 
that needs to be addressed. To this end, we present a very efficient linear 
inversion algorithm for local helioseismology. It is based on a subtractive op- 
timally localized averaging (SOLA) scheme in the Fourier domain, utilizing the 
horizontal-translation invariance of the sensitivity kernels. In Fourier space the 
problem decouples into many small problems, one for each horizontal wave vec- 
tor. This multi-channel SOLA method is demonstrated for an example problem 
in time-distance helioseismology that is small enough to be solved both in real 
and Fourier space. We find that both approaches are successful in solving the 
inverse problem. However, the multi-channel SOLA algorithm is much faster and 
can easily be parallelized. 
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1. Introduction 



The inverse problem of local helioseismology is to use measurements (e.g., wave 
travel-time shifts) to infer the physical conditions in the solar interior. For 
a recent review of local helioseismology, see, e.g., |Gizon, Birch, and Spruit| 
(2010). Inversions have been used to study flows and wave-speed perturbations 



et al. 



et al. 



2001 



around sunspots (e.g. Kosovichev 1996 Gizon, Duvall, and Larsen 
|Zhao, Kosovichev, and Duvall 2001[ Couvidat et al. 



20091, flows associated with the supergranulation (e.g. Kosovichev and 



2000 Jensen 



2004| |Gizon 



Duvall 19971 Woodard 2007 Jackiewicz, Gizon, and Birch 2008a), and global 



scale flows (e.g. Basu, Antia, and Tripathy 1999 Haber et al. 2002 Zhao and 



Kosovichev 2004 Gonzalez Hernandez et al. 2008) 



Essentially all inversions that have been employed in local helioseismology 
are linear. These inversions are based on the assumption of a linear relation- 
ship between perturbations to a reference model for the solar interior and the 
corresponding changes in the helioseismic measurements. The assumption of 
linearity is reasonable for inversions in the quiet Sun (e.g. Jackiewicz et al. 



2007b Couvidat and Birch 2009). Within the context of linear inversions, the 



two main approaches are optimally localized averages (OLA: Backus and Gilbert 



1968) and regularized least squares (RLS: Paige and Saunders 1982). In OLA 
methods, the goal is to produce spatially localized estimates of conditions in 
the solar interior while also controlling the associated random noise. In RLS 
methods, the goal is to produce a model of the solar interior that provides the 
best fit to the data under particular smoothness conditions. 

The first three-dimensional (3D) inversions in local helioseismology were based 
on the RLS formalism (iKosovichevl 119961 iCouvidat et al.\ 120051) and carried out 



using the LSQR algorithm (an iterative method, Paige and Saunders 1982 ). RLS 



corresponds to Tikhonov regularization (Tikhonov 1963) in the mathematical 



literature. This approach continues to be used extensively in time-distance he- 



lioseismology (e.g. Zhao, Kosovichev, and Duvall 


2001 Zhao and Kosovichev 


2004 Zhao et al. 2007). Jacobsen et al. 


(1999 


1 introduced the Multi-Channel 



Deconvolution (MCD) approach to solving the RLS equations. In MCD, the 
(assumed) horizontal-translation invariance of the kernels is exploited to decou- 
ple the full RLS problem into a set of easily solvable problems, one for each 
horizontal wave vector. This method has been used by, e.g., |Jensen, Jacobsen, | 



and Christensen-Dalsgaard 


(1998 


); 


Jensen et al. 


(2001 


); 


Couvidat et al. (2004) 



The 3D-OLA approach is computationally impractical for typical time-dis- 
tance inversions due to the size of the matrices involved, as we will show in 
Section [3] An improved OLA variant, termed Subtractive OLA (SOLA) and 
introduced by Pijpers and Thompson (1992), allows one to perform fewer ma- 



trix inverse computations, yet does not reduce the sizes of the matrices. SOLA 
corresponds to what is known as the method of approximate inverse in the 



mathematical literature on regularization of inverse problems (Louis and Maass 



1990 Schuster 2007). In problems where the kernel functions are separable as 
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products of functions of horizontal position and functions of depth, it is possible 
to reduce the 3D-SOLA problem to a set of 2D-SOLA problems followed by ID 



depth inversions (Jackiewicz et al. 2007a Jackiewicz, Gizon, and Birch 2008a). 



However, this separation is not always possible. 

In this article we will show an efficient Fourier method for carrying out 
3D-SOLA problems for local helioseismology. This method requires horizontal- 
translation invariance of the reference model for the solar interior and is based 
on an MCD approach. We will show that, unlike direct solution of the SOLA 
equations, this method is computationally feasible. We will also demonstrate that 
the multichannel approach is many orders faster than the real-space method with 
example computations. 



2. Setup of the Inverse Problem 



The solar interior is filled with numerous scatterers such as flows, magnetic fields, 
hot and cold spots, density and pressure anomalies, etc. The scattering mech- 
anism for each of these perturbations is physically different. Consider P such 
perturbations acting on the wave field. Not accounting for the magnetic field, 
the thermodynamic and flow perturbations would sum to P = 5 independent 
quantities. For example, the three components of flow velocity, sound speed, 
and first adiabatic exponent. 



In time-distance helioseismology (Duvall et al. 1993), the measurements con- 



sist of travel times between points and concentric annuli or between points and 
quadrants. These travel-time measurements are performed for different choices of 



Fourier filters (e.g. Gizon and Birch 2005 Jackiewicz, Gizon, and Birch 2008b). 



Taken together, we assume that we have M such measurements, each denoted by 
the running index a, where 1 < a < M. In this article we are concerned with the 
local helioseismology of a small patch of the Sun near disk center, and we make 
the approximation that sphericity can be ignored. Thus we adopt a Cartesian 
coordinate system: 



(r,z) = (x,y,z), 



(1) 



where r is the horizontal position vector on the solar surface and z is height. 
The statement that a number of scatterers acts on the wave field to create 



small shifts in travel times [5r a (r)] may be expressed as follows (e.g. Gizon and 
Bi7c^[2002l): 



5r a (r) 



d 2 r'dz K p{ r ' - r , z)8qp(r', z) + n a (r), 

13=1 



(2) 



where 5qp(r, z) represents the P perturbations in the various physical quantities 
that describe the solar interior, indexed by (3. The sensitivity of a travel-time 
measurement [<5r a (r)] to a localised change [5qp(r',z)\ is given by the travel- 
time sensitivity kernel [Kp(r' — r, z)\. For point-to-annulus or point-to-quadrant 
measurements [8T a (r)\ the position vector r usually denotes the center of the an- 
nulus (although there is some freedom in this convention). Note that in Equation 
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([2| we have explicitly assumed that the background solar model as well as the 
sensitivity kernels are invariant under horizontal translations. Sensitivity kernels 
result from a forward modeling under the single-scattering Born approximation 



( Gizon and Birch 2002 ) . The integral is taken over the volume of the Sun. The 



term n a (r) is the stochastic noise of the travel-time measurement [<5r a (r)] due 
to the forcing of waves by turbulent convection. The travel-time noise covariance 
matrix [Al has elements 



Kb{r t -Tj) = Cov [n a (r i ),n b (r j 



(3) 



Details about the computation of A can be found in Gizon and Birch ( 2004 ) 



The general OLA inversion problem for time- distance hclioscismology seeks 
to find an estimate of 8q a at any chosen target position [ro;^o]j given a set of 
St, K, and A. In other words, we are looking for a linear combination of the 
travel times such that 



JV M 



(4) 



= 1 a=l 



is an estimate of Sq a (rQ] z Q ). The weights w"(r.j — r ; z ) are the unknowns of the 
problem. In the above equation, N = (2n + l) 2 is the total number of horizontal 
position vectors [r$]. Throughout the article we assume that the travel times are 
given on a uniform square grid with sampling h x in both horizontal directions 
and with n x — n y — 2n + 1 pixels on each side. 



3. Subtractive Optimally Localised Averaging 

Using the Equations ^ and Q we have: 



r p 

6q™(r ;z ) = / dVd*^ 

® /3=1 



JV M 



^^<(r,-r„; Zo )^(r'-r„ Z ) 



N M 



|^»:(r,-r ;z K(r,). 
We can rewrite Equation ^ as: 



Sqr(r ;z ) = / d 2 r'dz1C«(r' -r ,z;z )5q a (r',z) 
J 

p 



5qp(r',z) 
(5) 

(6a) 



+ / dVd* V K1(r' -r ,z;z )6q p (r',z) (6b) 

JO a -i Q—L _ 



N M 



(6c) 



= 1 a=l 
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where the functions Kfp are averaging kernels given by 

N M 

/C^(r,z;z )=EE<( r - z o)^( r - r - z )' VjflGfl.P], (7) 



i— 1 a—1 



where /3 is a running index between 1 and P that labels the physical quantities. 
The term ( 6a ) on the right-hand side of Equation ^ is what we are searching 



for, i.e. the quantity Sq a convolved with the averaging kernel /C". If the averaging 
kernel /QJj(r, z; z ) is well localized in the horizontal (around r = and vertical 
(around z — zq) directions, we will recover a smoothed estimate of 5q a . 



The term (6b I is the leakage from the other perturbations /3 ^ a into Sq 1 ™ 
Ideally, one would like all averaging kernels JCg with /3 7^ a to be zero. 



The term (6c) represents the propagation of random noise from the travel 
times into the inverted Sq™ v . 

The SOLA method consists of searching for the inversion weights w"(r;zo) 
so that the averaging kernels ICp resemble user-supplied target functions [Tg Q ], 
while keeping error magnification as small as desired. This can be achieved by 
minimising the cost function 

r p 

X a (w a ;fi) = / d 3 ^ [iC a p(x;z )-T p a (x;z )] 2 (8a) 

® 13=1 
N M N M 

+ ^^2^2^2^2<(rt;z Q )A ab (r i -r J )w^{r j -z ). (8b) 

i=l a=l j=l 6=1 



Equation ([8j) has two main components: The first term (8a I, the "misfit", is a 
measure of how the averaging kernels [/Cg] match the target functions [Tg*]. To 
minimise the cross-talk, the target has only one non- vanishing component [T"], 
which we write as the product of a 2D Gaussian in the horizontal coordinates 
by a ID function of the vertical coordinate. Thus we write 



7?(x;zo)=Cexpl-^-)f(z-,Zo)6 Pa , V/3e[l,F], (9) 

where 5/3 a is the Kronecker 5- function. Typically, the function / peaks at a 
desired target depth z = z . The constant C is taken so that the spatial integral 
of Ta is unity. The parameter s controls the width of the Gaussian. 



The second term (8b I of Equation (|8j) is proportional to the variance [a„] of 



the random noise in g™ v due to the propagation travel-time noise: 

N M N M 

a l = HHlHl<( r ^ z o)Aab(r. i - r J )w^(r J ;z ). (10) 

i=l a=l j = l 6=1 

The trade-off parameter /i in Equation ^ is chosen to provide a satisfactory 
trade-off between the misfit and the noise; this choice is somewhat subjective. 
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The minimization is also subject to the constraints 

[ JC%(x;z )d 3 x = 6p a , V/3e[l,P], 
J 

which ensure that the inverted quantity is normalised appropriately. 



(11) 



4. Linear System of Equations 

The solution to the SOLA problem defined in Section [3] has been traditionally 
solved in real space for the one- and two-dimensional cases {e.g. Pijpers and 



Thompson 1994). Below we write the problem for the 3D case. 



We recast the optimization problem of Equation ([8]) subject to the constraints 
(111 by minimising the function 



(12) 



with respect to the vector of weights w a and a vector of Lagrange multipliers 
A. There are M x N unknown weights w"(rj; zq) and P unknown Lagrange 
multipliers A* 3 . The minimization consists of solving the following set of M x 
N + P equations: 







dw%(ri;z ] 



{C a } = 0, V(a,i) G [1,M] x [1,N], 



and 



ox? 



{C a } = 0, V/3e[l,P]. 



(13) 



(14) 



Equations ( 13 1 and ( 14 1 imply 



N M P 
3 = 1 6=1 P=l 

V(o,i) € [1,M] X [1,N], (15) 



and 



jV M 
3=1 6=1 



V/3G[1,P], 



(16) 



where we define 



r p 

® p=i 



+ /iA ab (r 4 - rj), 



(17) 
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C a p = [ K a (r-r l: z)d 2 rdz = [ K%{x)d?x, (18) 
J J® 

t2(r f ,z ) = [ KZ(r- ri ,z)T«(x;z )d 3 x. (19) 
Jo 



Equations (15) and (16) form a system of M x TV + P linear equations. The 



weights w"(ri] Zq) can be obtained by matrix inversion, or some equivalent linear 



solver. Inspecting Equations (17), (18) and (19), one sees that the dimension of 
the matrix that is to be inverted is (MN + P) x (MN + P). An advantage of the 
SOLA method (see Section [3]) compared to the OLA method is that the systems 
of equations for different target functions differ only in the right-hand sides and 



not in the matrix. Therefore the large matrix in Equation (17) has to be set up 



and inverted only once. For example, after inversion, we can infer the physical 
quantity at any depth using a new z ). 

Let us now discuss the computational cost of the SOLA scheme for typical 
time-distance inversions. The computational costs depends very much on the 
size of the matrix A to be inverted. Typically, A may have over 10 14 elements 
(over a petabyte!), which do not fit in computer memory. In practice, we do not 
solve the full problem, but truncate the sums over j in Equations ( 15 ) and ( 16 1. 
We consider a restricted number of convolution "shifts" [rj — (xj,Vj)] such that 



^shifts h x < Xj < n s hif tB h x , 
-n B hHtsh x < Vj < n s hihsh x , 



(20) 
(21) 



where h x is the horizontal sampling. The total number of shifts [iV s hitts = 
(2n s hift s + l) 2 ] must be much less than N so that the problem can now be 
solved. The minimum number of shifts that is acceptable depends on the size of 
the target function and the horizontal extent of the sensitivity kernels. In the 
example of Section [6j we have n s hifts = 45. The smaller the number of shifts, the 
worse the approximation of the problem and the worse the localisation and the 
noise of the answer. 

If one takes the same number of parameters as in the (2+l)D flow inversion 



of Jackiewicz, Gizon, and Birch (2008a) then P = 3 (the three components of 
velocity), n s hift s = 10, so that N sh [ {ts = 441, and M = 3 x 5 x 20 for three 
geometries (waves propagating in "outward minus inward" , "West - East" , and 
"North -South" directions), five ridges (/, pi, P2, _P3, and p4 modes), and 
twenty different radii. The problem would then require inverting a matrix of 
size w 100 000 x 100 000 that occupies tens of gigabytes of memory. Further- 
more, it is well-known that matrix inverse operations scale as 0(J\f 3 ), where 
Af = (NA1 + P) is the length of the matrix on one side. Even without memory 
issues, this calculation becomes intractable very quickly. 

Another issue that should be mentioned is the computation of the kernel- 



overlap integral in Equation (17). This computation is extremely expensive in 



real space. However it can be sped up very significantly by transforming to 
horizontal Fourier space, where it becomes a simple multiplication. Convolution 
operations are 0(N 2 ) in the real space but reduce to 0(Af\nAf) when performed 
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in Fourier space. In order to describe this step explicitly, let us define the dis- 
crete Fourier transform. Any function / and its Fourier transform / are related 
according to 

/>) = (||a£/(»-)e- lfc - r . (22) 

f(r) = *!£/(fc)e*- r , (23) 
fc 

where — 2n / (n x h x ) is the sampling in Fourier space. The horizontal wave 
vector fc takes the discrete values k pq — ph^e x + qh^i-y, where p and q are 
integers in the range [— n, n]. Using this definition of the Fourier transform, we 
obtain 

N P 

h l E E - ^z)K\{r x - n,z) = 

1=1 /3=1 

P 

(2nh k ) 2 £ e *-( r '- r ^ £ K a p *(k, z)K*(k, z), (24) 
fc p=i 

where we used the fact that the Fourier transform of K%{—r) is K%*(k) since 
Kg is real. In this form the kernel-overlap integral is computed much faster 



than in the Equation (17). The SOLA inversion examples presented later were 
computed using these equations. In order to avoid the edge effects resulting 
from the implicit periodicity assumed by the Fourier transform, we padded the 
sensitivity kernels and noise covariance matrices with zeros over a zone as wide 
as the size of the widest sensitivity kernel. 



5. Solution in Fourier Space 

In this section we fully exploit the horizontal-translation invariance of the sensi- 
tivity kernels and rewrite the entire problem in Fourier space. Using the definition 
of Equation (22), the Fourier transforms of Equations (15) and (16) are 

M P 

htNj2 A ab(k)w%(k ] z ) + 5 ktO Y,C a ^ = hl%(k-,z ), Va,fc (25) 

6=1 0=1 



and 



htN^C bf) vjZ{0;zo) = 5f Sa , 



V/3, 



(26) 



6=1 



where fc takes the discrete values k pq — ph k e x + qhke y , with p and q in the 
range [— n, n\. This set of equations can be written conveniently in matrix form 
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a.s 



if fc^O, 
if k = 0, 



h 2 NA(k)w a (k;z Q ) =t (k;z ), 



KNA(0) C 
C T 



w a (0;z ) 



h 2 k t a (0;z ) 
U a /(h 2 N 2 ) 



(27) 



where the vector t° (k) = [tf(k;zo) ^(fc^o) •■■ tfi(k', zq)] and the matrix 
A(fc) = A a b(k) have the following elements 



A ab (k) = (2tt) 2 / ^Kf(k,z)k^{k,z)dz + nK ab (k), (28) 

J ^bot 8=1 



K a a *{k,z)f«{k,z-z )&z, 



(29) 



where Zbot and z t op are the bottom and top heights of the computation box. 
Furthermore, U a = [6i a 5 2l a ••• 3p,a] T , w a (k) — [w1(k; z ) w^ik; z ) . . . 



w"j(k; zq)} , and C = [C a p] has elements given by Equation (18) 



In Fourier space the problem decouples into many small problems, one for each 
horizontal wave vector k. These small problems are completely independent and 
therefore can be solved in a parallel fashion. The solution w"(k) is constructed 
for each k separately. By analogy to the RLS Multi-Channel Deconvolution 



(Jacobsen et al. 19991, we call the current approach multi-channel SOLA or 



MCD SOLA. 

For each wave vector, the matrix to be inverted is much smaller than in the 
real-space case. For each wave vector k ^ 0, the matrix is of size M 2 . Taking 
the same parameters as in Section [4j the Fourier approach would only need 
441 inversions of matrices of size 300 x 300. This would result in an increased 
speed by more that five orders of magnitude over the real-space method for this 
realistic example. Note that there is no need to truncate the problem anymore 
("shifts = n). 

We provide below expressions for the averaging kernel [Equation Q] and the 
variance of the noise [Equation ( 10 1] in terms of the Fourier transform of the 
weights: 



M 

JC$(r, z; z ) = htN^2e ik ' r £ <(fe; zo)K$(k, , 

k a=l 



(30) 



M M 



a=l 6=1 fc 



(31) 



We emphasize that the averaging kernel is now computed as a matrix multipli- 
cation instead of a convolution. 
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Sensitivity Kernel Target Fimcti. 




x [Mm] x [Mm] 

Figure 1. Inputs to the example inversion. The point-to-point Born sensitivity kernel for 
sound speed is shown in the left column, and a target function with a full width at half 
maximum 20 Mm in the right column. The kernel in the x — y plane (top left) is the ID spatial 
integral of the kernel over depth. The black and white circles denote the two observation points, 
separated by 15 Mm. The depth slice in the lower-left panel is taken along the y = line. 
The horizontal cut of the target function (top right) is at a depth of 1 Mm. The depth slice 
(lower rig ht) is also along y = 0. The depth profile of the target was computed according to 
Equation \33\. 



The inferred solar property <5<z" lv at position (r, Zq) is 

M 

Sqr(r;z ) = NhiY / ^ r T,<*^ z o) Sfa ( k ) > ( 32 ) 
k a=1 

where <5f a (fc) is the Fourier transform of the travel-time maps. 



6. Example Inversion for Sound Speed 

6.1. Setup 

We now show a rather simple example of a time-distance helioseismic inversion 
to demonstrate the Multichannel SOLA method and compare it to its real- 
space counterpart. For simplicity, we will only consider one mean (mn) point-to- 
point travel-time measurement with the distance between the two observation 
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Figure 2. Model noise covariance matrix for mean point-to-point travel times with A = 15 
Mm for the example inversion. The left panel shows the noise covariance matrix in units of s 2 
as a function of horizontal coordinates. The right panel is a cut through the matrix along the 
y = line. The averaging time for this noise estimation is eight hours. 



points fixed at A = 15 Mm. For consistency with the measurement, we have 
computed a point-to-point Born-approximation sensitivity kernel according to 



Birch, Kosovichev, and Duvall (2004). This kernel gives the sensitivity of mean 
travel times to the sound speed perturbation [5c 2 /c 2 ]. No prior filtering has 
been done, i.e. the whole model power spectrum is used. We are not interested 
in specific types of kernels for this example problem, only a comparison between 
inversion methods and to prove that the MCD inversion works. This sound- 
speed kernel, denoted K^ n according to the conventions in Section [2j is shown 
in Figure [T] This kernel has 91 x 91 elements in the horizontal direction and 80 
elements in the vertical direction. 

The second input quantity kept fixed for our example inversion is the tar- 
get function, shown alongside the kernel in Figure [T] This 3D function has a 
Gaussian horizontal structure with a full width at half maximum of 20 Mm 
[sec Equation Q]. Since we are only working with one single sensitivity kernel 
(one A) in this example, it would be futile to attempt to obtain an averaging 
kernel peaked at some chosen z = z , since the kernel itself possesses no such 
depth properties. Therefore, again for simplicity, we choose a depth profile of the 
target function by horizontally integrating the sensitivity kernel at each depth 
coordinate to obtain a one-dimensional curve according to 

f(z;z ) = J d 2 rK^(r,z). (33) 

The ID-curve f(z; zq) is then combined with the horizontal Gaussian to con- 
struct the 3D target as in Equation |9]). This choice for /(z; z ) keeps the example 
as simple as possible. Note that the target in Figure [I] has a weak negative lobe 
beneath a depth of 10 Mm as a result of the depth profile of the sensitivity 
kernel. 

The final input quantity to the inversion is the noise covariance matrix defined 
in Equation ^ and denoted in this case as A mnjmn . We compute the covariance 



from the model power spectrum according to Gizon and Birch ( 2004 1 and show 
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the results in Figure [2j What this matrix tells us is how two mean point-to- 
point travel-time measurements are spatially correlated due to noise as the pairs 
of observation points are moved around with respect to each other. Note that 
for this case there is a significant correlation only when the measurements are 
made within about 10 Mm of each other. 

6.2. Comparison of the Real-Space and Fourier-Space Solutions 

We compute a set of real-space inversions and one Fourier inversion using the 
input quantities. The Fourier inversion is not computed in parallel for this ex- 



ample. For each inversion we generate a trade-off, or "L" curve (Hansen 1998) 



by choosing ten values of the parameter /j, [see Equation ( 28 )] . The values of \i 
are chosen to span the space of misfit and noise. One trade-off curve is generated 
for the Fourier inversion, but several are generated for the real-space inversion, 
each corresponding to a different number of shifts employed. The possible n s hifts 
range from 1 to 45, with 45 being the maximum due to the size of the kernel for 
this example (where n x — n y = 91). Since the MCD-SOLA inversion, in some 
sense, utilizes all possible shifts once, the real-space method with 45 shifts and 
the multi-channel method should agree. 

In Figure[3]we show the results for these inversions. The top panel of Figure|3] 
shows the trade-off curves, with red lines indicating the real-space inversion for 
varying shifts indicated by the numbers at the bottom of the curves. The thick 
blue line is the MCD inversion trade-off curve. These curves are typically plotted 
as the square of the random noise level versus the misfit on a logarithmic scale. 
We see that for increasing n s hift s the real-space inversion solution tends to the 
MCD solution. In fact, the L-curve for the 45-shift inversion falls on top of the 
one for the MCD inversion. Also shown in Figure [3] is a particular inversion 
weights f° r each inversion, chosen from the first (topmost) point on each 
trade-off curve when \x = 0.01 Mm -3 . These are the points where the averaging 
kernel and target match best, i.e. smallest misfit. This is a reasonable choice 
since our main concern here is not the noise, which is still quite small anyway. 
The last two weights in Figure [3] demonstrate what the L-curves already suggest: 
the solutions of the two types of inversions are perfectly comparable when we 
take the maximum allowable number of shifts in the real-space method. The 
weights for inversions with a smaller number of shifts are quite ill-behaved due 
to edge effects, and in practice one actually never uses all possible shifts since 
it is computationally impractical to do so. This suggests that standard 3D-OLA 
inversions might have undesirable properties in the solution due to the necessary 
truncation of the problem. Cuts through all weights are shown in Figure [4j 
reinforcing this point when only a subset of shifts is used. 

We recorded the computation time for each inversion. For the 45-shift case, 
the convolution matrix size is 8281 x 8281. The real-space inversion took two 
orders of magnitude longer to compute than the Fourier inversion (100 seconds 
compared to 1 second). This distinction only gets larger as the problem gets 
larger. Simply stated, the Fourier inversion takes a fraction of the time for 
small problems; for large problems, the real-space inversion is computationally 
intractable. 
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Figure 3. Example inversion using both the real-space and MCD methods and one single 
sensitivity kernel. The top panel shows trade-off curves where the noise unit is the fractional 
difference in sound speed [8c 2 /c 2 ]. The thin red lines are for the real-space inversion using an 
increasing number of shifts from top to bottom indicated by the numbers. The thick blue line is 
the Fourier-space inversion. The compute time for the 45-shift inversion was about a factor of 
100 larger than for the MCD inversion. The bottom panel corresponds to the inversion weights 
at trade-off parameter fi = 0.01 Mm -3 . This is the topmost point of each curve. Note that the 
size of the weights in the real-space SOLA inversion depends on the number of shifts, which 
is why the spatial scale changes, although the color scale is fixed throughout. 




-40 -30 -20 -10 10 20 30 40 

a: [Mm] 

Figure 4. Cut through each inversion weight in Figure [3] along the line y = 0. The curve for 
the 45-shift real-space inversion and the MCD-SOLA inversion are almost indistinguishable. 



SOLA: ms.tex; 14 September 2011; 0:35; p. 13 



J. Jackiewicz et al. 




x [Mm] x [Mm] 



Figure 5. Comparison of averaging kernel (left column) and target function (right column) 
from the example MCD-SOLA inversion. The target is the same one shown in Figure [I] The 
top panels are slices through the averaging and target functions at a depth of 1 Mm. The 
bottom panels are slices with depth through the averaging and target functions along the 
y = line. 

To show that the inversion does indeed work, in Figures [5] and [6] we provide 
comparisons between averaging kernel and target from the MCD solution. In the 
horizontal direction the agreement is quite acceptable, especially considering we 
have used only one input sensitivity kernel. Since the vertical profile of the input 
target function was constructed to match that of the sensitivity kernel, the good 
agreement there is expected. 

7. Discussion and Conclusions 

The example considered here is a very simple, toy inversion to demonstrate 
the usefulness of this new Fourier-based MCD-SOLA method. We have also 
experimented with a larger problem whereby we consider point-to-point mea- 
surements of various orientations of the observation points with respect to the 
x-axis. We input the same kernel as the one shown in this work, as well as 
horizontal rotations of it to match the measurements. Using the MCD method, 
we found that only five rotations, spaced evenly between and 90 degrees and 
keeping A fixed, are needed to find a very good averaging kernel that does 
not change with the addition of more rotations. The computation time for this 
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Figure 6. One-dimensional cuts through the averaging kernel and target function from Fig- 
ure[5]for our example MCD inversion. The left panel is a horizontal slice through each function 
at a depth of 1 Mm. The right panel is a plot along the x = y = line with depth. 



inversion was about 1.5 seconds. Had we attempted to solve the same problem 
with the real-space OLA inversion, in addition to consuming 40 gigabytes of 
memory, it would have taken weeks to compute. 

We also solved our toy problem with several kernels of various distances [A] . 
This allows one to obtain some resolving power in depth since the sensitivities 
differ. A target function was chosen as a 3D Gaussian peaked at a depth of 4 Mm 
beneath the surface. The MCD-SOLA inversion was able to find, as expected, an 
almost identical averaging kernel as the standard SOLA method. For a realistic 
application of the MCD-SOLA method with various targ et depths from 5 Mm 



Svanda et al. 



(2011), 



to the surface, we refer the reader to the recent work of 
who inverted for vector flows using synthetic travel-time observations as input. 
In conclusion, in this article we have extended what was originally done for 



RLS inversions (Jacobsen et al. 1999) to a SOLA inversion. A toy example 



inversion problem was solved with this new approach to compare and contrast 
to the more standard real-space SOLA method. The example proved that the 
MCD-SOLA works completely satisfactorily while the real-space counterpart 
may be intractable for all but the smallest problems. In fact, we demonstrated 
that for a realistic helioseismic problem, the MCD-SOLA method can be orders 
of magnitude more computationally efficient than the corresponding real-space 
method. We focused here on applications to time-distance helioseismology, but 
this approach is completely generalizable to any local helioseismic method re- 
quiring inversions, such as ring diagram analysis and acoustic holography. With 
the vast amounts of seismic data from the Solar Dynamics Observatory (SDO), 
it is imperative to have efficient and consistent local helioseismic OLA inversion 
procedures for studying the solar interior. 
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